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SECTION  I 


INTRODUCTION 


The  guidance  systems  for  modern  airborne  missiles  are 
increasing  in  complexity  with  the  introduction  of  onboard 
computational  power  in  the  form  of  minicomputers  and 
microcomputers.  The  traditional  form  of  guidance  systems 
for  airborne  missiles  employs  some  form  of  proportional 
navigation  (pro-nav),  a  guidance  scheme  developed  in  the 
^SO's  for  use  with  analog  equipment.  Modern  technology 
with  microminiaturization  of  digital  components  allows  a 
redefinition  of  the  guidance  schemes  to  utilize  the 
computational  power  of  current  state-of-the-art  equipment. 
Adaptive  control  systems  with  variable  gains  and  onboard 
processing  capabilities  provide  the  potential  for  increased 
performance  and  versatility  for  these  weapons.  The  price 
paid  for  the  increase  in  potential  is  the  development  and 
application  of  modern  optimal  control  and  estimation 
techniques  to  the  missile  guidance  problem. 

The  sensor  complement  employed  in  air-to-air  missiles 
is  usually  composed  of  rate  gyros  and  accelerometers  for 
determining  information  about  the  states  of  the  missile  and 
a  seeker  for  determining  information  about  the  target.  The 
seeker  may  be  classified  as  passive  if  it  provides 
measurements  of  line-of-sight  data  only,  and  active  or 
semi-active  if  it  provides  measurements  of  line-of-sight 
data  and  range  data.  Range  data  and  line-of-sight  data  may 
be  obtained  from  onboard  radar  while  line-of-sight  data  only 
may  be  obtained  from  an  infrared  sensor.  Many  short  range 
air-to-air  missiles  employ  passive  seekers  and  that  is  the 
type  of  seeker  to  be  employed  for  this  study. 

The  missile  model  to  be  employed  in  this  study  is  that 
of  a  generic  short-range  bank-to-turn  (BTT)  missile.  This 
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type  of  missile  employs  short  aerodynamic  surfaces  (wings) 
which  can  be  used  to  generate  large  aerodynamic  forces  in 
the  pitch  plane  of  the  missile.  For  steering,  the  missile 
rolls  to  orient  the  pitch  plane  in  the  direction  of  the 
desired  control  force  and  uses  the  lift  force  as  the  control 
force.  This  type  of  steering  differentiates  the 
bank-to-turn  missile  from  the  skid-to-turn  (STT)  missile 
which  does  not  employ  wings  to  produce  the  turning  maneuver. 
The  skid-to-turn  missile  does  not  need  to  orient  any  plane 
during  the  turning  maneuver,  but  it  also  cannot  generate  the 
turning  maneuver  forces  comparable  to  the  bank-to-turn 
missile. 

Proportional  navigation  control  laws  are  developed 
around  passive  seekers  where  the  objective  is  to  null  the 
rotational  rate  of  the  1 ine-of-s  ight  vector  from  the  missile 
to  the  target.  For  non-maneuvering  targets,  this  will  lead 
to  an  optimal  interception.  For  maneuvering  targets, 
pro-nav  control  laws  are  quickly  seen  to  be  suboptimal.  The 
use  of  modern  optimal  control  theory  is  ideally  suited  to 
this  problem.  In  the  majority  of  cases,  the  development  of 
the  control  lav  follows  the  form  of  the  linear  quadratic 
optimal  control  problem  in  which  the  control  is  assumed  to  be 
composed  of  a  linear  combination  of  the  states  and  the 
optimal  control  is  chosen  to  minimize  a  quadratic  function 
known  as  the  performance  index.  The  assumption  that  all  of 
the  states  are  measured  and  available  is  not  true  in  general 
and  a  computational  technique  must  be  developed  to  obtain 
the  desired  state  information  from  the  available 
measurements . 

An  estimation  algorithm  is  employed  to  convert  the 
available  measurement  information  into  the  desired  state 
information.  For  problems  governed  by  linear  state 
differential  equations  with  linear  relationships  between 
state  and  measurement  variables,  the  optimal  estimator  is 
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the  Kalman  filter  algorithm.  For  problems  with  non-linear 
state  differential  equations,  or  with  non-linear 
state-measurement  relationships,  or  hot'  the  extended 
Kalman  filter  is  the  most  widely  used  estimation  algorithm. 
While  these  are  powerful  estimation  algorithms  and  do 
perform  well  in  the  presence  of  noise  and  uncertainty,  most 
people  use  the  "cookbook  approach"  of  applying  the 
algorithms  without  understanding  the  theory  or  even  the 
results. 

The  Kalman  filter  employs  a  model  of  the  state 
propagation  process  to  determine  estimated  state  values  at  a 
point  in  time.  With  these  estimated  state  values,  an 
estimate  of  the  measurements  can  be  made  and  compared  to  the 
true  measurements  made  at  that  time.  The  discrepancy 
between  the  estimated  measurements  and  the  actual 
measurements  may  then  be  used  to  update  the  state  estimates 
to  a  more  correct  value. 

One  of  the  problems  with  the  application  of  the  Kalman 
filter  algorithm  is  that  the  characteristics  of  the  error 
associated  with  the  process  of  propagating  the  state 
estimate  (process  error)  and  the  error  associated  with  the 
taking  of  the  measurements  (measurement  error)  must  be 
specified  for  use  by  the  algorithm.  The  specification  of 
these  characteristics  for  a  certain  problem  configuration  is 
known  as  "tuning"  of  the  filter  and  may  lead  to  degraded 
performance  of  the  filter  under  another  problem 
configuration.  A  method  of  alleviating  this  tuning 
difficulty  is  to  allow  the  filter  to  be  self-tuning.  This 
process  is  known  as  adaptive  filtering  and  will  be  the 
subject  of  the  adaptive  filtering  section  of  this  report. 

In  summary,  the  application  of  modern  optimal  control 
techniques  and  the  attendant  estimation  techniques  to  the 
guidance  of  an  airborne  missile  pursuing  a  maneuvering 
target  provides  an  excellent  means  of  demonstrating  the 


application  of  modern  guidance  schemes.  In  addition,  the 
problem  provides  an  opportunity  for  a  comparison  of  modern 
guidance  schemes  to  traditional  guidance  schemes  and,  most 
of  all,  an  excellent  problem  for  development  and  testing  of 
new  guidance  and  estimation  techniques.  The  current  study 
places  perspective  upon  the  techniques  employed  in  modern 
missiles  while  providing  means  of  increasing  the  accuracy 
and  performance  of  the  missile  guidance  systems. 

It  should  be  emphasized  that  the  guidance  and  control 
algorithms  of  this  study  along  with  the  attendant  estimation 
algorithms  are  being  designed  with  the  idea  of  microcomputer 
or  minicomputer  implementation  onboard  the  missile.  For 
this  reason,  algorithms  of  the  recursive  form  are  emphasized 
rather  than  algorithms  which  require  large  memory  storage 
areas  or  extended  computational  time  per  iteration. 
Computational  efficiency  and  minimal  storage  requirements 
are  desirable  characteristics  of  the  optimal  guidance  scheme 
which  is  to  be  physically  realizable  onboard  the  missile. 


SECTION  II 


PROBLEM  DESCRIPTION 


This  section  of  the  report  provides  a  physical 
description  of  the  problem  to  be  studied  and  defines  the 
coordinate  systems  which  will  be  used  in  the  problem 
description.  Conventional  right-handed  coordinate  systems 
are  employed  for  the  problem  description. 

1.  Coordinate  Systems 

At  the  instant  the  airborne  missile  is  launched  from 
its  carrier  aircraft,  an  inertial  coordinate  system  is 
established  with  its  origin  at  the  center  of  mass  of  the 
missile,  with  x-axis  and  y-axis  in  the  horizontal  plane  and 
z-axis  vertical.  The  positive  x-axis  is  in  the  original 
direction  of  travel  for  the  missile,  the  z-axis  is  positive 
downward,  and  the  y-axis  completes  the  right-handed  system. 
This  is  the  base  coordinate  system  for  all  conversions. 

As  the  missile  moves  away  from  the  point  of  launch,  a 
set  of  missile  body  axes  moves  with  the  missile.  This  set 
of  axes  has  an  origin  always  at  the  center  of  mass  of  the 
missile,  the  x-axis  is  the  missile  longitudinal  axis  with 
the  positive  direction  forward,  the  y-axis  is  in  the  plane 
of  the  missile  wings  with  the  positive  direction  out  the 
right  wing,  and  the  z-axis  completes  the  right-handed 
system.  The  missile  body  axes  coincide  with  the  inertial 
axes  at  the  instant  of  launch,  but  then  the  body  axes 
translate  and  rotate  with  the  missile  whereas  the  inertial 
axes  remain  fixed.  The  missile  position  and  orientation  are 
described  in  terms  of  the  relationship  between  the  inertial 
axes  and  the  missile  body  axes. 

In  a  manner  similar  to  the  missile  body  axes,  a  set  of 
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target  body  axes  is  established  with  origin  at  the  center  of 
mass  of  the  target,  the  x-axis  is  the  longitudinal  axis  of 
the  target  with  the  positive  direction  forward,  the  y-axis 
is  in  the  plane  of  the  wings  with  the  positive  direction  out 
the  right  wing  of  the  target,  and  the  z-axis  completes  the 
right-handed  system.  The  target  body  axes  directions  rarely 
coincide  with  the  missile  inertial  axes  directions 
established  at  the  instant  of  launch.  The  target 
orientation  is  therefore  described  in  terms  of  the 
relationship  between  the  inertial  axes  and  the  target  body 
axes. 

Figure  1  provides  pictorial  representation  of  the  above 
information  for  a  number  of  air-to-air  engagements. 


Figure  1.  Engagement  Geometry  Definitions 


2. 


Missile  Mode  1 


The  missile  model  used  in  the  simulation  represents  a 
generic  bank-to-turn  missile.  This  is  a  highly 
maneuverable,  short  range,  air-to-air  missile  which  employs 
short  aerodynamic  surfaces  (wings)  to  produce  large 
aerodynamic  forces  in  a  plane  perpendicular  to  the  wings  to 
effect  a  turn.  This  means  that  a  turning  maneuver  by  the 


missile  requires  banking  the  missile  to  orient  the  plane  of 
the  aerodynamic  force  and  then  generating  the  aerodynamic 
force  to  cause  the  turn.  This  leads  to  the  common 
designator  of  bank-to-turn  for  this  category  of  missiles. 
This  type  of  missile  is  capable  of  generating  large  turning 
accelerations  for  maneuvers  (in  excess  of  100  Gs  at  moderate 
angles  of  attack),  and  care  must  be  exercised  or  the 
structural  limits  of  the  missile  will  be  exceeded.  In  this 
phase  of  the  study,  an  autopilot  furnished  by  the  USAF 
Armament  Laboratory,  Eglin  AFB,  Florida,  was  used  to 
transform  the  commands  from  the  guidance  system  to  the 
appropriate  bank-to-turn  maneuver  for  the  missile.  The 
missile  is  launched  from  the  carrier  aircraft  and  guidance 
commands  are  ignored  for  the  first  half-second  to  allow  the 
missile  to  clear  the  carrier  aircraft.  Guidance  commands 
are  then  obeyed  as  dictated  by  the  guidance  system.  The 
fuel  on  board  the  missile  lasts  for  approximately  two  and 
one-half  seconds. 

Following  the  thrusting  phase,  the  missile  follows  a 
ballistic  trajectory  utilizing  the  bank-to-turn  capability 
to  attain  the  target.  Sensors  on  board  the  missile  consist 
of  rate  gyros  and  accelerometers  to  detect  the  linear  and 
angular  accelerations  of  the  missile,  but  these  sensors 
provide  noisy  information. 

Additionally,  the  missile  has  a  sensing  system  for 
detecting  the  target.  A  passive  sensor  system  consisting  of 


an  infrared  seeker  is  used  to  provide  angular  position  data 
for  the  target  relative  to  the  missile.  Angular  rate  data 
was  not  used  for  this  study.  This  passive  seeker  does  not 
provide  information  directly  usable  by  the  guidance  system, 
so  an  estimation  process  must  be  employed  to  extract  the 
desired  information.  An  active  seeker  (one  which  provides 
range  data  in  addition  to  angular  data)  would  make  the 
estimation  process  much  more  accurate,  but  the  study  was 
restricted  to  use  of  the  passive  seeker,  as  the  use  of  tec 
active  seeker  was  prohibitive  in  both  weight  and  cost. 

3.  Target  Model 

The  target  used  in  this  report  is  a  "smart"  target 
which  employs  an  evasive  maneuver  in  an  attempt  to  escape 
the  incoming  airborne  missile.  The  target  maintains  a 
constant  speed  and  heading  until  the  missile  is  within  6000 
feet.  At  this  time,  the  target  begins  a  9-G  acceleration 
maneuver  in  a  direction  dictated  by  the  engagement  geometry. 
This  maneuver  is  continued  until  the  missile  is  within  1000 
feet.  At  this  time,  the  target  executes  a  "last  chance" 
maneuver  consisting  of  a  9-G  acceleration  vertically 
downward.  This  evasive  maneuver  has  been  tested  by  research 
at  Eglin  AFB  and  determined  to  be  the  type  of  maneuver  \;hich 
exercises  the  estimation  and  guidance  algorithms  in  the 
maximum  for  the  engagements  studied. 

Target  aerodynamics  are  not  modeled  in  this  study  as 
that  would  restrict  the  validity  of  the  study  to  the  type  of 
target  used  in  the  model.  Instead,  the  target  is  modeled  as 
a  rigid  body  whose  original  orientation  is  specified  by  the 
target  body  axes  orientation.  The  target  velocity  vector  is 
along  the  target  body  x-axis,  and  the  evasive  maneuver 
acceleration  is  along  the  target  body  z-axis.  Due  to  the 
short  time  periods  involved  in  these  air-to-air  combat 
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engagements,  the  target  thrust  acceleration  may  be  ignored, 
and  the  evasive  acceleration  is  obtained  from  aerodynamic 
forces  acting  upon  the  target.  The  9-G  evasive  acceleration 
was  chosen  as  a  typical  structural  limit  and  may  be  changed 
in  the  simulation  if  so  desired. 
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SECTION  III 


MATHEMATICAL  DEVELOPMENT 


This  section  of  the  report  describes  the  development 
of  the  differential  equations  of  motion  utilized  in  the 
problem  development.  In  addition,  the  equations  relating 
the  measured  quantities  obtained  from  the  sensors  to  the 
quantities  needed  by  the  guidance  system  are  provided. 

1.  State  Equation 

The  state  vector  differential  equation  of  motion  for 
an  airborne  missile  pursuing  an  airborne  target  can  be 
expressed  in  a  number  of  different  ways.  One  possible  means 
of  expression  is  in  terms  of  the  target  inertial  state 
vector  consisting  of  target  position,  velocity,  and 

acceleration  quantities,  and  the  missile  inertial  state 
vector  consisting  of  missile  position,  velocity,  and 

acceleration  quantities.  This  provides  a  total  of  18 
quantities  to  be  determined.  Choosing  the  missile 
acceleration  vector  as  the  control  to  be  specified  reduces 
the  number  to  15  quantities  to  be  determined.  The  guidance 
system  requires  full  knowledge  of  the  state  vectors,  so  the 
estimation  process  would  have  to  furnish  15  quantities  from 
its  input  sensor  information  (2  angular  measurements).  In 
order  to  reduce  the  dimensionality  of  the  estimation 
process,  note  that  the  missile  guidance  system  is  trying  to 
reduce  the  distance  between  the  target  and  the  missile  to  a 
minimum.  This  does  not  require  the  inertial  position  of  the 
missile  and  the  target,  but  the  relative  position  of  the 
target  to  the  missile.  Also,  the  inertial  velocity  of  the 
target  and  missile  are  not  required,  but  the  velocity  of  the 
target  relative  to  the  missile  is  required.  Lastly,  the 
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inertial  acceleration  of  the  target  and  the  missile  could  be 
replaced  by  the  relative  acceleration  of  the  target  to  the 
missile.  This  last  replacement  is  not  needed  as  the 
inertial  acceleration  of  the  missile  is  the  control  quantity 
which  will  be  specified  to  produce  the  state  quantities. 
Replacing  the  target  inertial  acceleration  with  the  target 
relative  acceleration  does  not  reduce  the  dimensionality  of 
the  problem.  Use  of  the  relative  position  and  velocity  of 
the  target  and  the  inertial  acceleration  of  the  target 
provides  a  problem  state  vector  with  nine  quantities  to  be 
determined  along  with  three  control  quantities  (the  missile 
inertial  acceleration)  to  be  specified. 

All  vectors  in  the  following  developments  are 
expressed  in  terms  of  components  along  the  inertial  axes 
which  were  established  at  the  instant  of  missile  launch. 
Let  R  denote  position  vector,  V  denote  velocity  vector,  A 

denote  acceleration  vector,  and  J  denote  the  jerk  vector 
(time  derivative  of  acceleration  vector).  Let  subscript  "T" 
denote  target  inertial  quantities,  subscript  "M"  denote 
missile  inertial  quantities,  and  subscript  "R"  denote  target 
quantities  relative  to  the  missile.  Once  these  symbols  have 
been  defined,  the  following  equations  governing  the  problem 
state  variables  may  be  written. 

R  =  V  -  V  =  V 
R  “T  ~M  R 

V  =*  A  -  A 
R  T  M 

A  =  J 
T  T 

In  order  to  continue  with  the  mathematical  description 
of  the  problem,  the  time  derivative  of  the  target  inertial 
acceleration  must  be  examined.  If  a  non-zero  value  is  used 
for  this  term,  then  these  equations  would  imply 


foreknowledge  of  the  target's  actions,  whereas  these  actions 
could  be  totally  random  and  not  describable  by  differential 
equations  at  all.  However,  an  assumption  can  be  made 
concerning  the  time  correlation  of  the  target  acceleration 
(a  random  variable).  The  success  of  this  practice  for  this 
particular  problem  did  not  warrant  using  an  assumed 
correlation  function  for  target  acceleration  with  the 
attendant  complexity.  Therefore,  the  time  derivative  of  the 
target  acceleration  should  be  set  equal  to  zero.  This 
action  does  not  preclude  target  acceleration  but  assumes 
that  the  target  acceleration  is  constant.  If  the  target 
acceleration  is  indeed  changing,  then  it  becomes  the  task  of 
the  estimation  process  to  adapt  to  the  new  acceleration  as 
rapidly  as  possible. 

Again,  using  the  inertial  axes  established  at  the 
instant  of  missile  launch,  define  the  inertial  components  of 
the  nine-element  state  vector  for  the  problem  in  the 
following  manner. 

XI  *  target  position  relative  to  missile  along  x-axis 

X2  *  target  velocity  relative  to  missile  along  x-axis 

X3  *  target  inertial  acceleration  along  x-axis 
Y1  *  target  position  relative  to  missile  along  y-axis 

Y2  *  target  velocity  relative  to  missile  along  y-axis 

Y3  *  target  inertial  acceleration  along  y-axis 
Z1  »  target  position  relative  to  missile  along  z-axis 

Z2  *  target  velocity  relative  to  m.ssile  along  z-axis 

Z3  *  target  inertial  acceleration  along  z-axis 

In  a  similar  way,  define  the  three-element  control 
vector  for  this  problem  in  the  following  manner. 

U1  »  missile  inertial  acceleration  along  x-axis 

U2  •  missile  inertial  acceleration  along  y-axis 

U3  *  missile  inertial  acceleration  along  z-axis 


With  the  state  vector  elements  and  the  control  vector 
elements  defined,  the  nine  scalar  differential  equations  of 
motion  for  this  problem  may  be  written  in  the  following 


manner  . 


XI 

= 

X2 

X2 

at 

X3 

-  U1 

X3 

= 

0 

Y 1 

3 

Y  2 

Y  2 

= 

Y3 

-  U  2 

Y3 

* 

0 

Z1 

= 

Z  2 

Z2 

S 

Z3- 

U  3 

Z3 

3 

0 

The  state  variable  differential  equations  of  motion 
may  now  be  written  in  the  standard  linear  equation  format  as 
follows . 


X  -  A  X  +  B  U 
Where : 

X'  -  [  X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3  ] 
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Note  the  natural  division  of  the  full  nine-state 
differential  equation  into  three  independent  three-state 
differential  equations. 


X  »  A1  X  +  B1  U1 


Y  *  A2  Y  +  B2  U2 

Z  -  A3  Z  +  B3  U3 
Where  : 

X'  =  [  X1.X2.X3  ] 

Y'  *  I  Y1.Y2.Y3  ] 

Z'  =  [  Z1.Z2.Z3  1 


A1  -  A2  =  A3  -  TO  1  0] 

lo  0  1 


The  state  variable  differential  equations  are  now  in 
the  standard  form  of  linear  differential  equations  with 
constant  coefficients  whether  the  standard  form  refers  to 
the  full  nine-state  system  or  to  one  of  the  three 
independent  three-state  systems.  The  state  differential 
equation  is  a  linear,  non-homogenous ,  vector  differential 
equation  with  constant  coefficients  and  standard  solution 
techniques  involving  the  state  transition  matrix  may  be 
applied.  Refer  to  Appendix  A  for  a  review  of  this  solution 
technique . 

2.  Measurement  Equations 

The  measurements  by  which  the  missile  control  system 
relates  to  the  real  world  consist  of  two  angular 
measurements  made  by  the  passive  seeker  on  board  the 
missile.  From  the  discussion  on  the  coordinate  systems  used 
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in  this  problem,  remember  that  an  inertial  coordinate  system 
is  established  at  the  instant  of  launch.  If  the  seeker  is 
an  inertial  platform  initialized  at  launch  with  measurements 
made  in  the  inertial  axes,  then  the  angles  are  exactly  the 
ones  needed  for  this  development.  If  the  seeker  is  a 
"strapdown  seeker"  where  the  measurements  are  made  in  body 
axes,  then  a  transformation  to  the  inertial  axes  is 
necessary.  The  following  development  assumes  that  the 
angular  measurements  are  made  in  the  inertial  axes. 

The  angular  measurements  are  made  by  an  infrared  "heat 
seeker"  on  board  the  missile  which  provides  the  angles 
defining  the  line-of-s ight  (LOS)  direction  from  the  missile 
to  the  target.  This  direction  is  specified  by  two 
angles  —  the  azimuth  angle  which  is  measured  in  the 
horizontal  plane,  and  the  elevation  angle  which  is  measured 
in  the  vertical  plane.  Figure  1  provides  a  pictorial 
definition  of  these  two  angles.  The  relationship  between 
the  state  variables  (relative  position,  relative  velocity, 
and  target  inertial  acceleration)  and  the  angular 
measurements  (inertial  azimuth  angle  and  inertial  elevation 
angle)  are  given  by  the  following  equations. 

Tan (  AZ  )  -  YI  /  XI 

2  2.5 

Tan(  EL  )  «  -  ZI  /  (  XI  +  YI  ) 

Where : 

AZ  is  the  current  LOS  inertial  azimuth  angle 
EL  is  the  current  LOS  inertial  elevation  angle 

XI  is  the  relative  position  along  the  inertial  x-axis 

YI  is  the  relative  position  along  the  inertial  y-axis 

ZI  is  the  relative  position  along  the  inertial  z-axis 

3.  Problem  Definition  Angles 

In  addition  to  the  measurement  angles  which  are  used 


in  the  estimation  algorithm  to  establish  the  relationship 
between  the  estimated  states  and  the  actual  states,  there 
exists  a  set  of  angles  which  are  used  to  identify  the 
particular  problem  geometry.  These  angles  are  the  angle  in 
the  horizontal  plane  defining  the  original  direction  of  the 
line-of-s ight  from  the  carrier  aircraft  to  the  target 
aircraft  (the  of f-bores  ight  angle),  and  the  angle  in  the 
horizontal  plane  defining  the  original  direction  of  travel 
for  the  target  aircraft  (the  target  aspect  angle).  As  the 
problem  is  extended  to  different  altitudes  for  the  initial 
positions  of  the  carrier  aircraft  and  the  target  aircraft, 
two  more  angles  will  become  necessary.  For  the  co-altitude 
engagement,  the  of f-bores ight  angle  and  the  target  aspect 
angle  are  sufficient  to  specify  the  engagement. 

These  angles  are  shown  on  Figure  1  where: 

TAZ  is  the  original  target  velocity  direction 

TOB  is  the  original  1 ine-of-s ight  direction 


SECTION  IV 


ESTIMATION  ERROR  ANALYSIS 


The  Kalman  filter  is  an  optimal  estimation  process  for 
a  totally  linear  system,  but  the  extended  Kalman  filter  and 
its  derivatives  are  suboptimal  when  applied  to  nonlinear 
systems.  If  the  propagation  equation  for  the  state  estimate 
is  known  exactly  ,  then  the  initial  uncertainty  indicated 
by  the  initial  value  of  the  state  estimation  error 
covariance  matrix  will  soon  decay  to  a  small  or  zero  value. 
If  the  state  propagation  equation  is  not  known  exactly,  then 
the  value  of  the  state  process  noise  covariance  matrix  will 
determine  if  the  state  estimation  error  covariance  will 
decrease  or  increase  as  the  state  is  propagated  forward  in 
t  ime  . 

The  common  tendency  on  all  of  the  estimation 
algorithms  tested  was  for  an  initial  increase  in  the  state 
estimation  error  covariance  during  the  period  that  the 
target  was  furtherest  from  the  missile,  followed  by  a  rapid 
convergence  toward  zero  as  the  missile  nears  the  target  and 
the  measurement  changes  become  large  in  comparison  to  the 
initial  measurement  changes. 

This  section  of  the  report  will  discuss  the  four  types 
of  estimation  algorithms  which  were  tested  and  the  five 
types  of  air-to-air  combat  engagements  which  were  studied  to 
analyze  the  performance  of  the  filters.  Results  will  be 
produced  which  will  show  the  time  histories  of  the  position 
error,  the  velocity  error,  and  the  acceleration  error  for 
each  engagement  and  for  each  filter. 

Table  1  identifies  the  individual  filters  by  number 
and  Table  2  identifies  the  engagements  by  number. 


FILTER  IDENTIFICATION  NUMBERS 


TABLE  1. 

FILTER  FILTER 

NUMBER  NAME 


1  EXTENDED  KALMAN  FILTER 

2  ADAPTIVE  KALMAN  FILTER 

3  ADAPTIVE  TIME  WEIGHTED  KALMAN  FILTER 

4  ITERATIVE  KALMAN  FILTER  (3  ITERATIONS) 


TABLE  2.  ENGAGEMENT  IDENTIFICATION  NUMBERS 


ENGAGEMENT 

RANGE 

ASPECT 

OFF-BORESIGHT 

( num) 

(ft) 

(deg) 

(deg) 

1 

13000 

90 

0 

2 

8000 

0 

0 

3 

11000 

90 

40 

4 

3000 

135 

0 

5 

3500 

135 

40 

1  . 


Filter  Comparisons 


This  section  of  the  report  will  present  a  comparison 
of  the  four  filters  for  each  engagement.  The  initial 
engagement  geometry  is  depicted  by  the  first  figure  of  each 
group,  followed  by  the  time  histories  of  the  position 
estimation  error,  the  velocity  estimation  error,  and  the 
acceleration  estimation  error.  In  each  case,  the  error  is 
calculated  by  subtracting  the  true  value  of  the  quantity 
(taken  from  the  simulation)  from  the  estimated  value  of  the 
quantity.  Note  that  the  position  refers  to  the  position  of 
the  target  relative  to  the  missile,  the  velocity  refers  to 
the  velocity  of  the  target  relative  to  the  missile,  and  the 


target.  A  final  comparison  of  the  filter  performance  in 
tabular  form  will  complete  the  analysis. 

a.  Engagement  1 

The  initial  problem  geometry  for  engagement  1  is 
presented  in  Figure  2.  The  carrier  aircraft  is  flying 
northward  at  a  speed  of  Mach  =  0.9  at  an  altitude  of  10,000 
feet  when  the  missile  is  launched.  The  target  aircraft  is 
located  directly  ahead  of  the  carrier  aircraft  at  the  same 
altitude  and  a  range  of  13,000  feet  at  launch,  but  is 
traveling  eastward  at  a  speed  of  Mach  «*  0.9. 

The  estimation  algorithm  is  initialized  with  the 
correct  state  values  as  these  are  available  from  the  carrier 
aircraft.  As  we  note  on  Figure  3,  the  position  error  is 
very  small  for  filter  1  until  the  missile  thrust  is  cut  off 
at  approximately  2.5  seconds.  After  this  time,  there  is  a 
moderate  buildup  of  position  estimation  error  until  the 

target  aircraft  begins  an  escape  maneuver  at  approximately 
3.5  seconds.  The  position  estimation  error  rapidly 

increases  due  to  the  inability  of  the  regular  extended 
Kalman  filter  to  track  rapidly  changing  parameters.  As  the 
position  estimation  error  progresses  to  a  peak  value,  the 
missile  is  drawing  nearer  to  the  target,  so  the  measurements 
are  becoming  more  usable  in  the  estimator.  This  is 

reflected  by  a  decrease  in  the  estimation  error  between  5 
and  6  seconds.  Just  after  6  seconds,  the  target  aircraft 

begins  its  final  maneuver  as  reflected  by  a  rise  in  the 
position  estimation  error.  The  missile  is  within  1000  feet 
of  the  target  at  this  time,  so  the  aircraft  maneuver  is  very 
influential  upon  the  estimation  algorithm  through  large 
changes  in  the  measurements.  The  estimation  algorithm 
reacts  to  this  good  data  by  rapidly  converging  toward  the 
true  position  with  the  estimated  position.  The  final  miss 


distance  is  5.23  feet  even  though  the  estimated  miss  is  much 
greater  than  this.  A  characteristic  of  all  of  these 
maneuvers  is  that  the  estimation  algorithm  will  produce  an 
acceptable  miss  distance  as  long  as  it  overestimates  the 
t ime-to-go . 

From  Figure  3,  the  velocity  estimation  error  closely 
follows  the  characteristics  of  the  position  estimation  error 
until  the  last  second  of  the  trajectory.  During  the  last 
second,  the  model  must  make  some  drastic  changes  in  order  to 
accommodate  the  measurements  it  is  receiving.  The  way  that 
the  position  error  is  decreased  is  by  a  large  change  in  the 
model  velocity.  This  leads  to  an  increase  in  the  velocity 
estimate  error  which  can  be  appreciable  near  the  end  of  the 
trajectory. 

From  Figure  4,  the  acceleration  estimation  error 
follows  the  velocity  estimation  error  characteristics  as  it 
should  since  the  position,  velocity,  and  acceleration  are 
related  through  the  model. 

Figure  5  presents  the  position  estimation  error  for 
filter  2,  an  adaptive  extended  Kalman  filter  which  uses  20 
measurements  to  produce  the  sample  measurement  noise 
characteristics.  One  important  item  to  note  is  that  the 
scale  of  the  position  estimation  error  plot  is  different  for 
this  filter  from  that  of  filter  1.  Due  to  the  smaller 
sample,  the  local  noise  characteristics  are  used  instead  of 
the  global  noise  characteristics.  This  allows  filter  2  to 
adjust  more  rapidly  to  changing  conditions  than  could  filter 
1.  This  fact  is  apparent  in  the  narrow  shape  of  the 
position  estimate  error  plot  indicating  rapid  adjustment  and 
in  the  lowered  peak  error.  In  a  similar  fashion,  the 
velocity  estimate  error  plot  and  the  acceleration  estimate 
error  plot  filter  2  have  reduced  values  from  those  of  filter 
1.  The  final  miss  distance  filter  2  is  4.47  feet,  not  a 
drastic  improvement,  but  the  performance  (as  measured  by 


miss  distance)  of  filter  1  was  pretty  good.  If  the  filter 
performance  all  along  the  trajectory  is  compared,  then 
filter  2  is  an  obvious  improvement  over  filter  1. 


The  use  of  local  noise  characteristics  instead  of  global  charac¬ 
teristics  provides  an  increase  in  overall  filter  performance  due  to 
the  ability  of  the  filter  to  react  rapidly  to  changing  conditions. 

Filter  3  was  created  especially  to  decrease  the  time  required  for  the 
filter  to  adjust  to  changing  conditions.  Filter  2  weighted  each  of 
the  measurements  equally  in  obtaining  noise  characteristics 
from  the  sample.  For  20  measurements  in  the  sample  window, 
this  means  a  weight  of  19/20  for  the  previous  measurement 
statistics  and  a  weight  of  1/20  for  the  latest  measurement. 
In  an  attempt  to  place  more  weight  on  the  latest 
measurement,  the  weighting  was  changed  to  4/5  for  the 
previous  measurement  statistics  and  1/5  for  the  current 
measurement.  This  represents  a  fourfold  increase  in  the 
importance  of  the  latest  measurement  over  that  of  filter  2. 

Figure  7,  the  position  estimate  error  plot  for  filter 
3,  is  a  very  convincing  argument  for  the  merit  of  this  idea. 
Notice  once  again  that  the  plot  scale  has  changed,  and  that 
the  plot  width  (indicating  response  time)  has  become  very 
narrow.  Peak  error  on  this  plot  was  about  one-fourth  of  the 
peak  error  of  the  corresponding  plot  for  filter  1  and  about 
one-third  of  the  peak  error  of  the  corresponding  plot  for 
filter  2.  Both  velocity  estimate  error  and  acceleration 
estimate  error  filter  3  are  smaller  than  those  of  filter  1 
or  2.  This  filter  does  such  a  good  job  all  along  the 
trajectory  that  it  is  no  surprise  that  the  final  miss 
distance  is  only  2.19  feet.  Note  that  if  miss  distance  is 
the  only  comparison  criterion,  then  all  of  the  estimation 
filters  provide  comparable  performance.  The  use  of  the 
error  time  history  allows  one  to  make  a  more  telling 
comparison.  Finally,  other  weighting  was  tried,  but  the 


weighting  given  above  provided  the  best  balance  of  rapid 
adjustment  with  small  overshoot. 

Aesthetically,  filter  4,  the  iterated  extended  Kalman 
filter,  is  very  pleasing  in  that  it  should  yield  performance 
improvement  on  each  iteration.  This  has  not  been  the  case 
in  that  the  filter  diverges  for  a  large  number  of 
iterations.  This  filter  is  still  undergoing  evaluation  but 
will  be  presented  here  for  comparison  purposes.  The  filter 
performance  and  final  miss  distance  are  comparable  to  those 
of  the  regular  extended  Kalman  filter  of  filter  1.  The 
smoothing  of  the  peaks  and  valleys  in  the  error  plots  points 
out  the  potentiality  of  this  technique.  More  work  is  being 
done  to  improve  the  performance  of  the  filter,  and  it  is 
hoped  that  this  filter  will  compare  favorably  to  filter  3  in 
the  end.  Figures  12,  13,  and  14  are  the  error  plots  for 
this  filter  for  engagement  1. 


Figure  6.  Position  Error  History  for  Engagement  1,  Filter  2 


Figure  7  . 
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The  initial  problem  geometry  for  engagement  2  is  given 
in  Figure  15.  The  carrier  aircraft  is  flying  northward  at 
an  altitude  of  10,000  feet  and  a  speed  of  Mach  =  0.9.  The 
target  aircraft  is  directly  ahead  of  the  carrier  aircraft  at 
the  same  altitude  and  at  a  range  of  3000  feet.  The  target 
aircraft  is  flying  northward  at  a  speed  of  Mach  =  C.9. 

The  performance  for  each  of  the  filters  is  similar  to 
that  of  engagement  1  with  the  notable  exception  of  filter  3. 
The  use  of  the  time-weighting  in  combination  with  the 
adaptive  feature  allows  this  filter  to  improve  its 
performance  such  that  the  maximum  position  error  for 
engagement  2  is  one-half  of  the  maximum  position  error  of 
engagement  1.  This  reduction  is  especially  notable  when  the 
performance  of  the  other  filters  is  similar  for  both 
engagements.  Figures  16  through  27  provide  the  time  history 
error  plots  for  position  estimate,  velocity  estimate,  and 
acceleration  estimate  for  this  engagement. 
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Figure  22.  Position  Error  History  for  Engagement  2,  Filt 
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Figure  25.  Position  Error  History  for  Engagement  2,  Filter  4 


Engagement  3 
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Figure  23  represents  the  initial  problem  geometry  for 
engagement  3.  The  carrier  aircraft  is  again  flying 
northward  at  10,000  feet  and  Mach  =  0.9  at  launch.  The 
target  aircraft  is  originally  located  on  a  line  40  degrees 
east  of  north  relative  to  the  carrier  aircraft  and  traveling 
in  an  eastward  direction  at  Mach  =  0.9.  The  initial  range 
is  11, 000  feet . 

Due  to  the  severe  starting  conditions,  the  missile 
guidance  system  and  the  estimation  algorithm  must  work  very 
hard  to  acquire  the  target.  Once  again,  the  results  of  the 
filter  performance  comparisons  are  similar  to  those  of 
engagement  1  with  the  regular  extended  Kalman  filter 
providing  the  worst  performance  followed  closely  by  the 
iterated  extended  Kalman  filter.  Next  comes  the  adaptive 
extended  Kalman  filter  with  approximately  half  the  peak 
error  of  the  other  two,  and  finally,  the  time  weighted 
adaptive  extended  Kalman  filter  with  a  peak  error  of  less 
than  one-sixth  the  peak  error  of  the  regular  extended  Kalman 
filter.  Figures  29  through  40  present  the  filter  error 
comparisons  for  engagement  3. 
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Figure  34.  Acceleration  Error  History  for  Engagement  3 
Filter  2 
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Figure  35.  Position  lirror  History  for  Engagement  3, 
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Figure  36.  Velocity  Error  History  for  Engagement  3,  Filter  3 
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Engagement  4 


Figure  41  represents  the  initial  problem  geometry  for 
engagement  4.  The  carrier  aircraft  is  again  traveling  north 
at  10,000  feet  and  a  speed  of  Mach  =  0.9  at  launch.  The 
target  aircraft  is  directly  ahead  of  the  carrier  at  the  same 
altitude  and  a  range  of  3,000  feet,  but  is  traveling  in  a 
southeastward  direction  at  a  speed  of  Mach  =  0.9.  This  is  a 
short  duration  missile  flight  due  to  the  initial  range. 
This  also  means  that  the  measurement  differences  are  large 
enough  to  cause  good  estimation  of  the  states.  All  of  the 
filters  performed  well  with  the  notable  fact  that  the 
maximum  position  estimate  error  of  the  four  filters  was 
experienced  by  filter  3.  This  maximum  error  was  less  than 
80  feet,  however,  and  filter  3  still  managed  to  attain  the 
smallest  miss  distance  in  the  end. 
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Figure  44.  Acceleration  Error  History  for  Engagement  4 
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Figure  47.  Acceleration  Error  History  for  Engagement  4 
Filter  2 
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e.  Engagement  5 

The  initial  problem  geometry  for  engagement  5  is 
presented  in  Figure  54.  Again,  the  carrier  aircraft  is 
flying  northward  at  10,000  feet  and  a  speed  of  Mach  =  0.9  at 
launch.  The  target  aircraft  is  also  at  10,000  feet  but  is 
located  on  a  line  40  degrees  east  of  north  relative  to  the 
carrier  aircraft.  The  target  aircraft  is  traveling 
southeastward  at  Mach  =  0.9.  The  initial  range  is  3500 
f  eet . 

This  is  another  short  range  problem  with  severe 
starting  conditions.  All  of  the  filters  performed 
acceptably  in  this  condition  which  might  cause  problems  with 
conventional  proportional  navigation  systems.  Surprisingly, 
the  t ime-weigh t e d  adaptive  extended  Kalman  filter  performed 
poorest  of  the  group  (but  still  performed  acceptably)  while 
the  iterated  extended  Kalman  filter  performed  best  of  the 
group . 
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SECTION  V 


CONCLUSIONS 


Four  forms  of  Che  extended  Kalman  filter  have  been 
investigated  for  five  different  launch  conditions  in  a 
simulated  air-to-air  combat  scenario.  All  of  the  estimation 
algorithms  performed  in  an  acceptable  manner  in  that  all 
filters  allowed  the  missile  to  come  within  ten  feet  of  the 
target  and  that  is  considered  a  hit.  The  comparison  of  the 
filter  performance  along  the  entire  trajectory  provides  more 
insight  into  the  filter  performance  than  comparing  miss 
distance  alone.  The  regular  extended  Kalman  filter 
consistently  performs  the  poorest  when  maximum  estimation 
error  and  rapidity  of  adjustment  to  changing  conditions  are 
the  judging  criteria.  The  use  of  an  adaptive  filter  which 
uses  a  local  window  of  measurements  to  compute  the 
measurement  noise  statistics  significantly  improves  the 
performance  of  the  regular  extended  Kalman  filter.  The  use 
of  time  weighting  of  the  measurement  data  again 
significantly  improves  the  filter  performance  in  reducing 
the  maximum  estimation  error  and  also  in  adaptation  to 
changing  conditions.  The  iterated  extended  Kalman  filter 
did  not  perform  as  well  as  expected,  but  additional  work  in 
this  filter  is  expected  to  be  profitable.  The  iterative 
nature  of  the  filter  is  expected  to  smooth  the  error  history 
and  provide  better  estimates.  At  the  current  time,  the  use 
of  the  extended  Kalman  filter  which  uses  adaptation  to  local 
conditions  with  time  weighting  to  improve  the  rapidity  of 
adjustment  to  changing  conditions  is  highly  recommended. 

One  note  of  caution  for  the  casual  reader  these  filters 
have  been  used  with  a  linear  quadratic  control  law  where  the 
determination  of  the  final  time  is  paramount  in  determining 
the  optimal  control.  The  best  filter  will  perform  poorly  if 


coupled  with  a  poor  control  system.  The  future  work  of  this 
group  in  exploring  the  options  for  the  control  system  hold 
great  promise  for  improving  the  overall  system  performance 
and  thereby  improving  the  estimation  process. 


APPENDIX  A 


DISCRETE  STATE  PROPAGATION  EQUATIONS 


Let  the  following  expression  represent  the  vector 
differential  equation  for  a  system  governed  by  a  set  of 
linear,  first  order,  ordinary  differential  equations. 

X  -  A  X  +  B  U 


In  the  equation  above,  the  coefficient  matrices  A  and 
B  may  be  constant  matrices,  matrices  which  are  functions  of 
time  explicitly  ,  or  matrices  which  are  functions  of  X  as 
well  as  time.  For  this  application,  the  .matrices  are 
constant  coefficient  matrices. 

The  control  vector  U  may  be  time  varying,  but  for  this 
application,  the  control  vector  will  be  considered  piecewise 
constant.  This  means  that  the  control  value  will  be  held 
constant  between  commands  and  that  the  control  reaches  the 
commanded  value  immediately.  This  is  an  obvious 
approximation,  but  this  is  the  development  of  a  model  for 
the  estimation  process,  so  the  approximation  is  acceptable. 

The  solution  of  the  state  propagation  equation  is  most 
easily  obtained  through  use  of  the  principle  of 
superposition. 

Let  the  state  at  some  time  t(k)  be  propagated  from  the 
state  at  some  time  t(k-l).  Let  X( k) *Y (k) +Z ( k)  where  Y(k) 
represents  the  propagation  of  the  initial  conditions  from 
t(k-l)  to  t(k)  with  zero  forcing  function  and  Z(k) 
represents  the  state  response  at  t(k)  to  the  forcing 
function  between  t(k-l)  and  t(k)  calculated  from  zero 
initial  conditions. 

First,  consider  the  Y  response. 


Y(k)  =  e 


A  (t(k)-t(k-l  )  ] 


Y  ( k-1  ) 


Let  F(k,k-1)  designate  the  state  transition  matrix 
which  is  represented  by  the  exponential  matrix  term. 

A  [ t ( k) -t ( k-1 )  ] 

F(k, k-1 )  =  e 

Now,  the  solution  to  the  linear,  homogenous,  ordinary 
vector  differential  equation  may  be  written  in  terms  of  the 
state  transition  matrix  in  the  following  manner. 

Y(k)  -  F(k, k-1 )  Y(k-l) 

Now,  consider  the  response  to  the  forcing  function 
independently  of  the  response  to  the  initial  conditions. 
Again,  let  Z(k)  denote  the  response  at  time  t(k)  due  to  any 
forcing  function  between  t(k-l)  and  t(k).  For 
simplification,  we  will  consider  the  forcing  function  B  U  to 
have  a  zero  value  everywhere  within  the  interval  t(k-l)  to 
t(k)  except  for  a  small  time  interval  "dt"  located  at  time 
t(j)  between  t(k-l)  and  t(k).  The  response  at  t(k)  will 
therefore  be  the  response  to  an  impulse  occurring  at  t(j). 

Z  -  A  Z  +  B  U 

Since  this  equation  begins  with  zero  initial  conditions,  the 
response  will  remain  zero  until  reaching  the  time  t(j)  when 
the  forcing  function  becomes  non-zero.  The  time  interval 
"dt"  over  which  there  is  a  non-zero  forcing  function  is 
considered  to  be  so  small  that  the  differential  equation 
governing  this  portion  of  the  total  response  may  be 
represented  by  the  following  equation. 


At  the  end  of  the  interval,  "dt",  the  forcing  functon 
is  removed  and  the  Z-response  at  this  time  is  given  by  the 
following  equation. 

Z  (t(j)  +  dt]  =  B  U  [ t  ( j )  ]  dt 

The  interval,  "dt",  over  which  the  impulsive  forcing 
function  is  defined  is  so  small  that  the  result  may  be 
interpreted  as  an  instantaneous  change  in  the  state  occurring 
at  time  t(j).  The  response  at  time  t(k)  may  now  be  obtained 
by  propagating  this  initial  condition  from  time  t(j)  to  time 
t(k)  with  no  forcing  function. 

Z ( k )  *  F ( k , j )  Z(j)  =  F(k, j)  B  U(j)  dt 

If  there  are  two  small  intervals  located  at  times  t(i) 
and  t(j),  over  which  there  is  an  impulsive  forcing  function, 
then  the  total  response  at  time  t(k)  is  the  sum  of  the 
individual  responses  at  time  t(k). 

Z ( k )  -  FCk.i)  B  U(i)  dt(i)  +  F(k,j)  B  U(j)  dt(j) 

For  many  small  intervals  of  non-zero  values  of  the 
forcing  function,  the  total  response  at  t(k)  may  be  written 
as  a  summation  of  the  individual  responses  to  the  individual 
impu 1 se  s  . 

Z(k)  -  T  F ( k , i )  B  0( i)  dt(i) 

i 

As  the  series  becomes  large  and  the  size  of  the 
intervals  becomes  small,  then  the  series  becomes  an 
integral.  In  this  way,  we  may  consider  any  continuous 
function  between  t(k-l)  and  t(k)  to  be  composed  of  an 


infinite  number  of  small  intervals  with  the  function 
continuous  within  each  small  interval  but  varying  from  small 
interval  to  small  interval. 


Z(k) 


Jt  F<k’ 


t)  B  U(t)  dt 


The  total  solution  at  time  t(k)  may  now  be  obtained  from 
superposition. 


X(k)  -  Y(k)  +  Z(k) 


X(k)  -  F(k, k-1 )  X(k-l)  + 


,t(k) 

F ( k , t )  B  U(t)  dt 
t(k-l) 


If  we  now  assume  instantaneous  response  of  the  control  and 
assume  that  the  control  is  held  constant  in  the  interval 
between  t(k-l)  and  t(k),  the  discrete  times  when  the  state 
is  available,  the  equation  above  may  be  simplified. 


X(k)  -  F(k,k-1)  X(k-l)  + 


,t(k) 

F(k, t )  B  U(k-l)  dt 
t (k-1 ) 


Let  us  investigate  the  quantity  in  brackets.  Assume 
that  the  state  propagation  process  occurs  in  steps  of 
constant  time  increments  given  by  DT  *  t(k)-t(k-l).  The 
limits  on  the  integration  can  be  made  to  be  "0"  and  "DT”  for 
this  process.  All  appearances  of  t(k)-t(k-l)  which  happened 
in  the  earlier  development  may  now  be  replaced  with  DT . 


F(k, k-1 )  -  e 


(A  DT) 


The  state  transition  matrix  F(k,k-1)  is  a  constant 
matrix  for  a  constant  coefficient  matrix  ,"A",  and  for  a 
constant  interval,  "DT".  This  constant  state  transition 


matrix  will  be  identified  with  a  single  letter  "F"  from  this 
time  on. 


Define  a  control  transition  matrix,  G(k,k-1),  which 
may  be  used  to  show  the  influence  of  the  constant  control  in 
the  interval  t(k-l)  to  t(k)  upon  the  state  response  at  time 
t  (k)  . 


G(k,k-1) 


t  (k) 

F(k , t )  B  dt 
t(k-l) 


Note  that  the  functional  form  of  the  state  transition 
matrix  "F"  is  used  in  the  evaluation  of  the  control 
transition  matrix  "G",  not  the  constant  "F".  Also  note  that 
for  a  constant  interval  "DT",  the  control  transition  matrix 
will  be  a  constant  matrix.  The  single  letter  "G"  will  be 
used  to  refer  to  the  constant  control  transition  matrix.  *3 

n 

The  values  of  F  and  G  may  be  precomputed  and  stored  as 
they  are  constant  matrices.  Thus,  the  continuous 
differential  equations  governing  the  state  propagation  may 
be  replaced  with  the  discrete  recursive  state  propagation 
equations . 

X  -  A  X  +  B  0 


j 


becomes 


X(k)  -  F  X(k-l)  +  G  O(k-l) 


I 

This  is  the  form  to  be  used  in  subsequent  ! 

developments. 


I 


I 


APPENDIX  B 


LINEAR  QUADRATIC  CONTROL  LAW 

This  section  of  the  report  considers  the  development 
of  the  control  law  which  provides  optimal  control  to  a 
system  trying  to  minimize  a  chosen  performance  index. 

The  control  law  will  be  developed  for  a  system  which 
obeys  a  set  of  linear,  first  order,  ordinary  differential 
equations  with  constant  coefficient  matrices.  Following  the 
notation  used  in  the  equations  of  Appendix  A,  the  linear 
differential  equations  of  motion  may  be  represented  as 
follows . 

X  -  A  X  +  B  U 

Beginning  at  some  time  T,  the  problem  is  to  determine 
the  control  law  which  is  used  to  specify  the  control  "U" 
such  that  some  measure  of  performance  C'J",  the  performance 
index)  may  be  minimized.  For  this  particular  problem,  the 
performance  index  is  chosen  as  a  combination  of  the  final 
state  and  the  integral  control  power  required  to  produce  the 
final  state. 


J  -  .5  XF'  S  XF  + 


/TF 
.5  D 


'(t)  W  U(t)  dt 


The  matrices  S  and  W  are  positive-definite,  weighting 
matrices  with  S  being  the  final  state  weighting  matrix,  and 
W  being  the  control  power  weighting  matrix.  Both  of  these 
matrices  are  constant  matrices  for  this  problem.  The  symbol 
XF  refers  to  the  state  value  at  the  time  TF. 

Now,  considering  that  at  some  time  T,  the  current 


state  is  known,  the  problem  is  to  determine  U  which  produces 
the  optimal  trajectory  based  upon  the  performance  index 
while  obeying  the  differential  equations  of  motion. 

An  augmented  performance  index  is  formed  by  appending 
the  differential  equation  constraints  to  the  original 
performance  index  through  the  use  of  a  costate  vector,  "P". 


H  »  .5  U'W  U  +  P'(A  X  +  B  U  -  X) 


AJ 


.5  XF' 


S  XF  + 


dt 


In  the  expressions  above,  "H",  the  augmented  integrand 
of  the  performance  index  is  known  as  the  numerical 
Hamiltonian.  The  augmented  performance  index,  "AJ",  can  now 
be  minimized  with  no  constraints,  and  this  is  equivalent  to 
minimizing  the  original  performance  index  with  the 
differential  constraints. 

Set  the  first  variation  of  the  augmented  performance 
index  equal  to  zero,  and  obtain  the  necessary  conditions  to 
produce  a  minimum  of  the  augmented  performance  index  and  a 
minimum  of  the  original  performance  index.  Note  that  the 
original  performance  index  was  a  quadratic  form  with 
positive  definite  weighting  matrices,  so  only  a  minimum 
exists  for  this  problem.  Determining  the  necessary 
conditions  will  therefore  completely  solve  the  problem. 

Setting  the  first  variation  of  the  augmented 
performance  index  equal  to  zero  produces  the  following 
conditions . 

Boundary  condition: 

(XF'  S  -  PF']  DXF  +  HF  DTF  -  0 
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State  differential  equations: 


X  -  A  X  +  B  U 

Costate  differential  equations: 

P  -  -  A'  P 

Optimality  condition: 

-1 

U  -  -  W  B'P 

The  boundary  condition  furnishes  a  final  value  of  the 
costate  variable  in  terms  of  the  final  value  of  the  state 
variable.  For  now,  the  final  time  will  be  considered  fixed 
and  the  allowable  change  in  the  final  time,  "DTF" ,  will  be 
zero.  If  the  optimality  condition  is  used  to  eliminate  the 
control  from  the  state  and  costate  differential  equations,  a 
classic  two  point  boundary  value  problem  will  result  with 
the  values  of  the  state  known  on  the  initial  time  boundary 
and  the  value  of  the  costate  known  on  the  final  time 
boundary.  Most  problems  of  this  nature  must  be  solved  by  an 
iterative  procedure,  but  this  particular  problem  possesses 
an  analytical  solution  in  terms  of  the  final  time. 

Refer  to  the  mathematical  development  section  of  this 
report  for  the  development  of  the  state  differential 
equations . 


ft 


The  definition  of  the  control  weighting  matrix  as  a 
diagonal  matrix  in  combination  vith  the  diagonal  form  of  the 
final  state  weighting  matrix  allows  one  to  express  the 
performance  index  for  the  full  nine  state  problem  as  the  sum 
of  three  independent  performance  indices,  one  for  each  of 
the  coordinate  axes.  This  follows  from  the  fact  that  the 
state  differential  equation  for  the  full  nine  state  problem 
could  be  expressed  in  three  independent  sets,  one  for  each 
of  the  coordinate  axes. 


R) 


8 


fs 


A: 

* 
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AAA 


A . 


m.  3  «  * 


•  ft. 


X'  -  [  XI , X2 , X3  ] 


X  -  A1  X  +  B1  01 


J1  -  .5  XF'  SI  XF 


-TF 

+  J  * 


5  Ol'(t)  w  Ul(t)  dt 


Following  the  usual  procedure  for  minimizing  a 
performance  index  subject  to  a  set  of  differential 
constraints,  we  will  form  an  augmented  performance  index  by 
introducing  the  costate  vector  "P"  which  will  have  the  same 
dimension  as  the  state  vector  "X".  The  augmented 
performance  index  "AJ1"  may  be  written  in  the  following 


manner . 


HI  -  .5  w  U1  +  P'(A1  X  +  B1  01  -  X) 


AJ1  -  .5  XF'  SI  XF  +  /  Hl(t)  dt 


The  quantity  "HI"  in  the  expression  above  is  the 

numerical  Hamiltonian.  Performing  the  perturbation 
mathematics  about  the  optimal  trajectory  leads  to  the 

following  necessary  conditions  for  the  optimal  path. 

X  »  41  X  +  B1  01 

P  -  -  Al'P 

01  -  -  l  Bl'P 
w 

PF  -  SI  XF 


Consider  the  analytical  integration  of  the  costate 


t  *  i 


**•  V 


differential  equations  which  are  linear,  homogenous 
differential  equations  with  constant  coefficients. 


P'  -  [  P1,P2,P3  ] 


P  -  -  Al'P 

PF  =* 

SI  XF 

PI  - 

0  : 

P1F 

P2  - 

-  PI  : 

P2F 

P3  - 

P2  : 

P3F 

XI F 

0 

0 


Integrating  the  differential  equations  from  the  known 
time  T  to  the  unknown  time  TF  yields  the  following 
expressions  for  the  costate  vector. 

PI  -  X1F 

P2  .  ■  X1F  (TF-T) 

2 

P3  -  X1F  1  (TF-T) 

7 

The  analytical  solutions  for  the  costate  differential 
equations  are  now  known  in  terms  of  the  parameters  of  final 
time  and  final  position.  While  these  parameters  have  not 
yet  been  evaluated,  continue  to  express  the  rest  of  the 
solution  in  terms  of  these  two  parameters.  Using  the 
optimality  condition,  the  optimal  control  can  be  expressed 
in  terms  of  the  parameters  as  follows. 
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becomes  : 


U1  -  1  X1F  (TF-T) 
w 


The  optimal  control  is  now  expressed  in  terms  of  the 
parameters  of  final  time  and  final  position.  The  use  of  the 
state  transition  matrix  approach  will  allow  further 
simplification  of  this  problem. 

/TF 

F(TF,t)  B  Ul(t)  dt 
T 

Where : 

01 ( t )  -  l  X1F  (TF-t) 
v 


Since  the  analytical  form  of  the  state  transition 
matrix  is  required  for  the  integration,  the  following 
evaluation  is  helpful. 


F(T , t ) 


A  (T-t) 
e 


yields  : 


F  -  A  F 


F ( t , t )  »  I,  the  identity  matrix 


Using  the  values  of  the  matrices  for  this  particular 
problem  yields  the  following  state  transition  matrix  for 
this  problem. 
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F ( TF , t ) 


(TF-t)  1  (TF-t) 

7 

1  (TF-t) 


The  analytical  expression  for  the  state  transition 
matrix  may  now  be  used  to  propagate  the  state  from  the 
current  time  "T"  to  the  final  time  "TF".  The  expression  for 
X1F  is  the  following. 

2  3 

XI F  -  XI  +X2  (TF-T)  +  X3  1  (TF-T)  -  X1F  1  (TF-T) 

7  Jw 

The  first  three  terms  in  the  equation  above  arise  from 
propagating  the  current  state  to  the  final  time  using  the 
state  propagation  matrix.  The  last  term  arises  from  the 
integration  of  the  optimal  control  influence  upon  the  state 
between  the  current  time  and  the  final  time.  Since  the 
difference  between  the  current  time  and  the  final  time  is 
the  true  parameter  instead  of  final  time,  define  the 
"time-to-go"  parameter,  TTG  ■  TF-T,  and  use  this  expression 
in  the  following  work.  Solving  for  the  final  state  from  the 
equation  above  yields  the  following  expression. 


XI F  -  Cl  [  XI  +  X2  TTG  +  X3  1  TTG  ] 

7 


Where 


Cl  -  (3  w)  /  (3  w  +  TTG  ) 


Therefore 


U1  ■  C2  [  XI  +  X2  TTG  +  X3  1  TTG  ] 


Where  : 


3 

C2  -  (3  TTG )  /  (3  w  +  TTG  ) 

Nov,  from  similarity,  the  optimal  control  expression 
for  the  other  two  directions  can  be  inferred  from  the  one 
dimensional  example.  The  optimal  control  for  each  of  the 
three  directions  may  be  expressed  as  follows 

2 

U1  -  GN  [  XI  +  X2  TTG  +  X3  1  TTG  ] 

7 

2 

U2  -  GN  [  Y1  +  Y2  TTG  ♦  Y3  1  TTG  ] 

7 

2 

03  -  GN  [  Z1  +  Z2  TTG  +  Z3  1  TTG  ] 

7 

Where : 

3 

GN  -  (3  TTG)  /  (3  w  +  TTG  ) 

The  optimal  control  for  this  problem  has  now  been 
expressed  in  terms  of  the  current  state  and  the  single 
parameter  "TTG",  the  time-to-go.  The  specification  of  TTG 
in  combination  with  the  current  state  knowledge  will  provide 
a  unique  value  for  the  optimal  control. 

The  determination  of  a  value  for  the  parameter  "TTG" 
is  a  study  in  itself.  Even  assuming  perfect  knowledge  of 
the  state  at  the  current  time,  how  does  one  predict  the  time 
required  for  the  airborne  missile  to  reach  a  maneuvering 
target  whose  actions  are  unpredictable?  The  solution  for 
TTG  is  usually  obtained  from  an  assumptiom  about  the  future 
target  acceleration  and  from  the  knowledge  of  the  missile 
acceleration  along  its  body  x-axis  which  is  uncontrollable. 
In  the  simplest  case,  assume  that  the  target  acceleration  is 


zero  and  that  the  missile  closing  velocity  is  constant.  The 
value  of  TTG  is  then  obtained  by  dividing  the  target 
distance  from  the  missile  by  the  closing  velocity.  Other 
means  of  obtaining  TTG  involve  approximating  the  closing 
acceleration  as  a  constant  value  and  solving  the  resulting 
quadratic  for  TTG.  More  exotic  means  of  calculating  TTG 
exist  but  all  of  the  methods  suffer  from  the  same  malady. 
The  value  of  TTG  obtained  is  valid  only  for  the  condition 
used  to  obtain  TTG.  Changes  in  the  thrust  acceleration  of 
the  missile  or  changes  in  the  target  acceleration  will  cause 
discontinuities  in  the  TTG  calculations.  For  this 
particular  problem,  major  discontinuities  occur  when  the 
missile  exhausts  its  fuel  supply,  and  when  the  target  begins 
its  evasive  maneuver.  A  variation  of  the  constant  closing 
acceleration  method  is  used  in  determining  TTG  for  this 
report,  but  this  technique  is  not  advocated  above  the 
others.  Future  studies  will  explore  additional  information 
concerning  the  determine tion  of  TTG  as  this  is  a  most 
important  parameter. 

With  the  determination  of  TTG,  the  optimal  control  law 
is  completely  specified  in  terms  of  the  current  state  vector 
information.  One  should  remember  the  assumptions  made  in 
the  determination  of  the  guidance  law  and  in  the 
determination  of  the  parameters  which  provide  a  unique  value 


APPENDIX  C 


DEVELOPMENT  OF  AN  ESTIMATION  ALGORITHM 


This  section  of  the  report  deals  with  the  development 
of  an  estimation  algorithm,  a  computational  technique  for 
extracting  desired  information  from  available  information. 
In  Appendix  B,  the  optimal  control  law  was  developed 
presuming  that  all  information  regarding  the  state  vector 
was  available.  The  estimation  process  is  used  to  furnish 
the  necessary  state  information  to  the  guidance  system.  The 
combination  of  a  Linear  Quadratic  control  law  with  an 
estimation  algorithm  which  assumes  Gaussian  noise 
characteristics  produces  a  control  system  known  as  a  Linear 
Quadratic  Gaussian  (LQG)  control  system. 

For  control  of  an  airborne  missile  with  a  passive 
seeker,  the  available  sensors  provide  information  about  the 
inertial  acceleration  of  the  missile  and  the  angles  defining 
the  1 ine-of-s ight  direction  from  the  missile  to  the  target. 
All  measurements  are  corrupted  by  some  form  of  additive 
white  noise,  and  there  is  uncertaintity  in  the  initial  state 
vector  quantities  at  the  beginning  of  the  engagement.  From 
the  noisy  measurements  and  the  uncertain  starting 
conditions,  the  estimation  algorithm  is  expected  to  furnish 
to  the  control  law  the  target's  position  and  velocity 
relative  to  the  missile,  and  the  inertial  acceleration  of 
the  target — the  full  nine  state  vector. 

1.  Extended  Kalman  Filter 

For  fully  linear  systems,  systems  obeying  linear  state 
propagation  equations  and  linear  measurement  equations,  the 
optimal  estimation  algorithm  is  the  Kalman  filter.  For 
those  systems  with  nonlinear  state  propagation  equations  or 
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with  nonlinear  measurements,  the  standard  estimation  process 
is  the  extended  Kalman  filter,  a  suboptimal  filter  which 
linearizes  the  nonlinear  part  of  the  process  about  a  "best 
estimate".  The  extended  Kalman  filter  is  the  base  algorithm 
for  this  report  and  all  other  estimation  algorithms  will  be 
compared  to  the  extended  Kalman  filter. 

The  estimation  algorithm  relies  upon  a  model  of  the 
process  by  which  the  state  vector  describing  the  system  is 
propagated  in  time.  This  problem  is  assumed  to  obey  the 
linear  model  developed  in  the  Mathematical  Development 
section  of  this  report  with  an  extra  term  to  account  for  any 
unmodeled  terms,  incorrectly  modeled  terms,  or  any  other 
errors  in  state  propagation  which  can  be  included  in  the 
category  of  "process  noise".  The  true  state  vector 
propagation  is  represented  by  the  linear  vector  differential 
equation  below. 

X-AX  +  BU  +  C 

Where  : 

X  is  an  N-state  vector 

A  is  an  NxN  matrix  of  constants 

B  is  an  NxM  matrix  of  constants 

0  is  an  M- control  vector 

C  is  an  N-vector  of  noise 

The  integration  of  the  state  differential  equations 
utilizing  the  state  transition  matrix  provides  the  following 
discrete  recursive  equation  governing  the  true  state 
propagation. 


m 
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X(k+1 )  -  F  X(k)  +  G  U(k)  +  q(k) 

Where : 

X ( k )  is  the  state  vector  at  the  k-th  time  point 
F  is  the  constant  state  transition  matrix 
G  is  the  constant  control  transition  matrix 
U(k)  is  the  control  vector  at  the  k-th  time  point 
q(k)  is  the  process  noise  for  the  k-th  interval 

At  this  time,  note  that  if  the  exact  initial  state 
vector  is  known,  if  the  exact  control  vector  is  known,  and 
if  the  exact  state  propagation  model  is  known,  then  this 
model  would  yield  the  true  state  vector  with  no  additional 
work.  The  only  inputs  required  by  the  model  would  be  the 
time  for  which  the  state  vector  is  required. 

The  fault  with  this  line  of  reasoning  is  that  there 
are  no  "exact"  quantities  in  the  real  world.  For  problem 
solution  in  the  real  world,  we  must  make  a  "best  guess"  or 
estimate  of  the  initial  state  vector,  produce  a  control 
vector  based  upon  the  estimated  state  vector,  and  propagate 
the  state  estimate  to  a  new  time  based  upon  an  "assumed 
form"  or  model  of  the  state  propagation  process.  At  the  new 
time,  information  in  the  form  of  real  world  measurements  is 
used  to  correct  or  "update"  the  state  estimate.  The 
estimation  process  consists  of  a  continuous  sequence  of 
propagation  and  updating  of  the  state  estimate. 

Let  XT  denote  the  true  value  of  the  state  vector,  let 
XH  denote  the  best  state  estimate  using  all  available 
information,  and  let  XB  denote  the  propagated  state  before 
including  new  measurement  information. 

XB(k+l )  -  F  XH(k)  +  G  U(k) 

XT(k+l )  -  F  XT(k)  +  G  U(k)  +  q(k) 

Remember  that  the  true  state  XT  is  not  known  and  that 
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the  XT  propagation  equation  above  is  written  for  comparison 
with  the  XB  propagation  equation  only. 

The  real  use  for  the  equations  above  is  the 
determination  of  a  state  estimate  error  propagation 
equation.  Define  the  error  terms  as  follows. 

XBe  -  XB  -  XT 

XHe  -  XH  -  XT 

Subtracting  the  true  state  propagation  equation  from 
the  estimated  state  propagation  equation  yields  the 
following  state  estimate  error  propagation  equation. 

XBe(k+l)  -  F  XHe(k)  -  q(k) 

The  error  propagation  equation  above  is  not  really 
used  to  propagate  the  state  estimate  error  as  we  would  use 
the  error  value  and  make  the  estimate  correct  if  possible. 
The  error  propagation  equations  are  useful  for  propagating 
error  bounds  which  can  be  established.  Next,  the  sign  of 
the  state  estimation  error  is  unknown,  so  the  sign  of  the 
error  bound  is  also  unknown.  The  use  of  a  covariance  matrix 
solves  some  of  these  problems.  Define  the  covariance  matrix 
for  an  error  vector  as  the  expected  value  of  the  outer 
product  of  the  error  vector  with  itself. 

PB  -  E{  XBe  XBe'  > 

PH  -  E {  XHe  XHe'  > 

In  the  expressions  above,  "PB"  is  the  covariance 
matrix  for  the  propagated  state  estimate  error,  and  "PH"  is 
the  covariance  matrix  for  the  updated  state  estimate  error. 
The  use  of  the  covariance  matrices  relieves  us  from  the 
necessity  of  guessing  the  appropriate  signs  of  the  error 
terms  as  the  diagonal  terms  of  the  covariance  matrix  are  the 
squares  of  the  error  terms. 


The  use  of  the  covariance  matrix  also  relieves  us  froo 
the  task  of  guessing  the  initial  estimation  errors.  Assume 
that  the  initial  estimation  errors  are  uncorrelated  with 
zero  expected  values.  This  allows  us  to  set  the  initial 
values  of  the  off-diagonal  terms  in  the  covariance  matrices 
to  zero  and  the  diagonal  terms  to  the  square  of  the  initial 
error  bounds.  Remember  that  the  state  estimation  error 
covariance  matrix  is  a  measure  of  the  confidence  which  can 
be  placed  in  the  state  estimate  as  indicated  by  the  growth 
or  decline  of  the  estimation  error  bounds.  The  equation  for 
the  propagation  of  the  state  estimation  error  covariance 
matrix  is  obtained  as  follows. 

PB  -  E {  XBe  XBe'  > 

Where : 

XBe ( k+1 )  -  F  XH(k)  -  q(k) 

Thus  : 

PB (k+ 1 )  -  F  PH(k)  F'  +  Q(k) 

Where : 

E{  F  XHe  XHe'  F'  }  -  F  PH  F' 

E {  F  XHe  q'  >  -  0 

E {  q  XHe'F'  }  -  0 

E  {  q  q  '  >  ■  Q 

Where  Q  is  the  process  error  covariance  matrix. 

The  initial  state  vector  can  be  estimated  and  the 
propagation  equation  used  to  propagate  the  estimate  to  a  new 
time.  At  the  same  time,  the  initial  state  estimation  error 
bounds  can  be  established  and  the  error  bounds  used  to 
initialize  the  state  estimation  error  covariance  matrix. 
The  estimation  error  covariance  matrix  may  be  propagated  to 
the  new  time  along  with  the  state  estimate.  One  requirement 


for  this  propagation  is  that  the  process  error  covariance 
matrix  "Q"  must  be  defined.  An  estimate  must  be  made  of  the 
error  between  the  model  equation  used  to  propagate  the  state 
estimates  and  the  real  world  propagation  of  the  true  states. 
This  error  may  be  used  to  establish  an  error  bound  for  the 
process  error  and  the  error  bound  may  then  be  used  to  create 
the  process  error  covariance  matrix.  For  this  problem,  the 
process  error  covariance  matrix  was  constant  throughout  the 
engagement.  Note  that  the  process  error  covariance  matrix 
always  increases  the  estimation  error  covariance  matrix. 

Once  the  estimate  of  the  state  vector  is  propagated  to 
a  new  time,  the  relationship  of  the  state  vector  to  the 
measurement  vector  can  be  used  to  estimate  what  the 
measurements  should  be.  Comparison  of  the  estimated 
measurements  to  the  true  measurements  provides  a  means  of 
correcting  or  updating  the  state  estimates.  Let  KT  denote 
the  true  or  actual  real  world  measurements  and  let  MB  denote 
the  estimate  of  the  measurements  based  upon  the  propagated 
state  estimate  and  the  state-measurement  relationship  model. 
Let  "MB  -  h(XB)"  be  the  model  of  the  state-measurement 
relationship  and  let  "MT  ■  h(XT)  +  r"  be  the  actual 
state-measurement  relationship.  The  quantity  "r"  is  the 
measurement  noise  term  similar  to  the  process  noise  term  in 
the  state  propagation  equations.  The  measurement  error  is 
defined  in  a  manner  similar  to  the  state  error. 

MBe  -  MB  -  MT 

MBe  *  h(XB)  -  h(XT)  -  r 

Linearize  the  last  equation  about  the  propagated  state 
estimate  where  XT  -  XB  -  XBe.  This  will  give  an  equation 
relating  the  measurement  residual  or  error  to  the  propagated 
state  estimate  error. 

MBe  -  HX  XBe  -  r 
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The  incorporation  of  the  measurement  data  into  the 
estimation  process  will  allow  us  to  correct  or  update  the 
state  estimate  in  order  to  make  the  measurement  estimate 
conform  more  closely  to  the  actual  measurement.  A 
relationship  must  be  established  between  the  measurement 
residual  and  the  updated  state  estimate  in  order  to 
accomplish  this.  For  obtaining  a  linear,  unbiased  state 
estimate,  the  updated  state  estimate  will  be  formed  from  a 
linear  combination  of  the  propagated  state  estimate  and  the 
measurement  residual.  The  updated  quantities  are  given  by 
the  following  equation. 
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Ff  i  if.,  i .  i .  pjr  r  rj  -v  ?.  ^  a  *rv.' 


XH  -  XB  -  KG  (MB  -  MT ) 

XHe  -  XBe  -  KG  MB  e 

*  (I  -  KG  HX)  XBe  -  KG  r 

PH  -  (I  -  KG  HX)  PB  (I  -  KG  HX)'  +  KG  R  KG' 

Note  the  introduction  of  the  gain  matrix  "KG"  which  is 
the  Kalman  gain  matrix.  The  matrix  "KG"  is  an  NxM  matrix. 
The  discussion  below  is  concerned  with  the  determination  of 
a  value  for  KG. 

The  manner  in  which  KG  is  determined  specifies  the 
estimation  technique.  If  KG  is  chosen  to  minimize  the 
estimation  error  in  some  way,  then  the  estimator  is  a 
minimum  error  estimator.  If  the  value  of  KG  is  chosen  to 
minimize  the  estimation  error  covariance  in  some  way,  then 
the  estimator  is  a  minimum  variance  estimator.  This  report 
is  concerned  with  the  minimum  variance  type  of  estimator. 

Take  a  first  variation  of  the  PH  equation  with  KG  as 
the  perturbation  variable.  Setting  the  first  variation 
equal  to  zero  will  yield  the  following  expression  for  KG. 

KG  -  PB  HX'  (HX  PB  HX'  +  R)  1 

Note  the  pleasing  result  — that  the  value  of  KG  can  be 
determined  from  the  propagated  state  estimate,  propagated 
state  estimate  error  covariance,  and  the  measurement  error 
covariance.  This  means  that  KG  can  be  determined  explicitly 
from  propagated  and  measured  quantities.  There  is  no 
iteration  necessary  in  determining  KG. 

Substitution  of  this  expression  for  KG  into  the 
equation  for  PH  allows  an  alternate  expression  for  PH  as 
given  below. 


PH  -  (I  -  KG  HX)  PB 


This  expression  for  PH  is  computationally  more 
efficient  than  the  previous  expression,  but  the  insight  that 
PH  is  a  symmetrical  positive  semidefinite  matrix  is  not 
available  from  the  latter  expression.  Likewise, 
periodically  one  should  perform  a  check  that  the  latter 
means  of  calculating  PH  does  retain  the  positive 
semidef initne s s  of  PH. 

The  following  provides  a  flowchart  for  the  extended 
Kalman  filter  algorithm  which  has  just  been  developed. 


Precompute : 


Initialize : 


Begin : 


Propagate 


Measure : 


Compute : 


Compute 


Update 


F  -  F ( DT , 0 ) 
G  -  G(DT, 0  ) 

XH ,  PH  ,  Q ,  R 


T(k+I)  -  T(k)  +  DT 
XB(k+l)  -  F  XH(k)  +  G  U(k) 
PB(k+i )  -  F  PH(k)  F'+  Q(k) 

MT(k+l ) 

MB(k+l)  -  h(XB(k+l ) ) 

MBe  -  MB  -  MT 
HX(XB) 

KG  -  PB  HX'  (HX  PB  HX'  +  R) 

XH(k+l)  -  XB ( k+ 1  )  -  KG  MBe 
PH  -  (I  -  KG  HX)  PB 


119 


2.  Adaptive  Kalman  Filter 

In  the  development  of  the  extended  Kalman  filter,  it 
was  noted  that  potential  difficulty  exists  in  the  manner  of 
specification  of  the  process  error  covariance  matrix  "Q", 
and  the  measurement  noise  covariance  matrix  "R".  Adaptivn 
estimation  algorithms  attempt  to  alleviate  the  problems  b\ 
allowing  the  filter  to  adjust  to  the  unknown  noise 
statistics.  Ideally,  adaptive  filtering  would  allow  the 
estimation  process  to  account  for  nonlinearities,  for  error 
sources  not  included  in  the  model,  and  for  all  forms  of 
noise.  Practically  speaking,  this  is  not  realizable,  but 
the  adaptation  should  improve  the  filter  operation  and  allow 
an  improvement  in  the  estimation  of  the  states. 

The  measurement  residuals  provide  an  excellent  source 
of  information  for  the  estimation  of  the  measurement  noise 
statistics.  A  moving  window  sample  of  measurement  residuals 
will  be  used  to  produce  an  unbiased  estimate  of  the 
measurement  noise  covariance  matrix.  The  following 
relationships  will  be  used  for  calculating  the 
characteristics  of  the  measurement  noise. 

Measure:  MT(j):  j«l,nm 

Calculate:  MB(j):  j*l,nm 

Determine:  MBe(j):  j*l,nm 

Define :  Cl*l /nm 

Define:  em  ■  Cl  £  MBe(j) 

j 

Define:  C2*l/(nm-l) 

Define:  CS  ■  C2  T*  [MBe(j)-em]  [MBe(j)-em]' 

**j 

Remember:  CM  -  HX  PB  HX'  +  R 

Approximate:  R  »  CS  -  HX  PB  HX' 

The  above  expressions  allow  the  online  determination 
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of  the  measurement  noise  statistics  based  upon  the 
measurement  residuals  and  the  assumption  that  the  sample 
covariance  can  be  used  to  approximate  the  global  covariance. 
The  validity  of  this  approximation  is  borne  out  in  the 
filter  comparisons  vhere  the  adaptive  filter  consistently 
improves  upon  the  performance  of  the  regular  extended  Kalman 
filter.  The  only  fault  that  one  could  find  with  the 
approximation  is  that  the  accuracy  is  dependent  upon  the 
size  of  the  sample  used  in  the  computations.  In  almost  all 
cases,  a  sample  window  of  20  measurements  provided  the  best 
performance,  so  the  adaptive  filter  window  size  was  fixed  at 
this  value. 
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